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Abstract. We study the crossover of a finite one-dimensional (ID) bosonic 
ensemble from weak to strong interactions in harmonic traps and multi-well 
potentials. Although these systems are very common in experimental setups 
and have been studied theoretically, an analytical description is lacking. We 
perform Diffusion Quantum Monte Carlo calculations which we show to be in good 
agreement with results from analytical functions that we construct to describe 
these systems. For the harmonic trap we use the correlated-pair wave function 
which we introduced in [l] considering here much larger atom numbers, going 
beyond the few-body ensembles studied in pQ. We also investigate double and 
triple wells, changing correspondingly the uncorrelated part of the Ansatz to 
describe efficiently the single-particle behaviour. On-site effects beyond mean-field 
and standard Bose-Hubbard calculations that appear in densities being captured 
by our analytical functions are explored. 
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1. Introduction 

One-dimensional many-body problems have a long tradition in theoretical physics 
since a substantial number of them allows for an analytical description mainly by 
employing the Bethe-Ansatz. A seminal example is the Lieb-Liniger gas [5] for 
bosons interacting via a contact potential of arbitrary strength in the absence of 
an external potential with periodic boundary conditions (and extendable also to 
hard- wall boundaries [3 ). The existence of an exact solution for such a general 
many-body problem has lead to further studies, also trying to extend this idea to 
the presence of external traps. One other important peculiarity of ID gases is the 
existence of the so-called Tonks-Girardeau gas [2], where hard-core bosons interacting 
with infinite strength can be mapped to identical non-interacting fermions, since the 
infinite repulsion is mimicking the Pauli-exclusion principle. 

In the last decade, the excitement about ID physics has received significant 
attention by corresponding experimental realizations, since the preparation of cold 
gases in quasi- ID traps with strong transversal confinement has become possible. The 
Tonks-Girardeau gas [5] and also its counterpart for the strongly attractive excited 
state, the super Tonks-Girardeau gas [6] have been observed, and numerous other 
ID phenomena like the pinning transition [7] with a weak optical lattice have been 
explored, extending also to nonequilibrium quantum dynamics and thermalization 
processes [HJ . All these would not be possible without the high degree of experimental 
control of the trapping geometry via external fields [5] as well as the interaction of 
ultracold atomic ensembles via Feshbach [10 or confinement induced resonances 
a tool specific to quasi-lD waveguide-like systems. The precision of measurements 
has been recently extended to single-cite addressing in optical lattices [12], while the 
deterministic preparation of the system has reached the possibility of loading an exact 
number of atoms into the trap [T31 [Tl] . 

Theoretical studies of ID systems have been mainly focusing on the uniform and 
the harmonic trap case. Especially Monte Carlo simulations in configuration space 
have been performed using trial functions of Jastrow type as a guiding for the diffusion 
|15j . The common Bijl- Jastrow trial function involves two variational parameters: one 
for the two-body interaction part and one for the single-particle part. On the other 
hand, one-dimensional lattices have been studied typically within the Bose-Hubbard 
model |16) . which is a tight-binding approach and in the usual form cannot capture 
on-site effects in the wave-function as long as unperturbed Wannier on-site orbitals 
are used. Numerical approaches which demonstrate localization and derealization 
mechanisms, as well as on-site features like fragmentation and multi-particle effective 
interactions have been developed |17j . 

In this work we perform a comparative study between quantum Monte-Carlo 
calculations and results from analytical wave-functions for the description of ground 
state energies and correlation properties of bosonic systems in a harmonic and multi- 
well traps. Our approach initiates from the exact solution of the two-body problem 
in the trap |18j and the generalization in terms of a correlated-pair wave-function 
(CPWF) which we have proposed in [TJ in order to cover the few-body (up to 6 atoms) 
case in the harmonic potential. We extend the calculations we have performed in [TJ 
the many-body regime (up to 50 atoms) in the trap performing a diffusion Monte Carlo 
(DMC) simulation which provides us with the numerically-exact energy of the system 
which we compare with the corresponding numerical integration of the CPWF. We 
also extend the CPWF Ansatz to cover multi-well potentials (double and triple well) . 
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For the construction of a suitable parameter-free Ansatz here we use localized Wannier 
functions of Gaussian profile for the single particle part of the wave function to describe 
the arrangement of the particles in the multi-well trap, and hypergeometric functions 
inspired from the two-body problem in the trap to include pair-correlation and on-site 
effects. A Variational Monte Carlo (VMC) approach allowing for one or two variational 
parameters has been also employed without substantial differences with respect to the 
obtained energy values compared to our parameter-free Ansatz. We demonstrate 
on-site features appearing in multi-well potentials for several observables (one-body 
density and pair- correlation functions) captured by our analytical functions. 

This article is organized as follows: in Section II we introduce our setup, its 
Hamiltonian, as well as our Ansatz. In Section III, we compare results for the energy 
from DMC simulations with these from the trial function CPWF. We also show one- 
and two-body density functions, focusing on on-site effects. Finally we summarize our 
results and give an outlook in Section IV. 

2. Setup and Ansatz 

2.1. Hamiltonian for harmonic and multi-well traps 

For the investigation in particular of ID systems one should take into account the 
experimental conditions and their impact on the collisional properties of the atoms. 
Experimentally, the standard method to create quasi-lD tubes is using a very strong 
laser field for the transversal directions compared to the lateral one [6] . This way 
the trap becomes highly anisotropic with the characteristic transversal length scale 

a±_ = yj -p^-j- being much smaller than the longitudinal one <j|| = yl A ^ where uj± 

is the transversal (longitudinal) harmonic confinement frequency. Then the transverse 
degrees of freedom are energetically frozen to their ground states and the effective ID 
interaction strength reads [TTj : 



where do is the 3D s-wave scattering length. 

We will study here two systems. First as a standard benchmark example a ID 
harmonic trap where the rescaled ID N-body Hamiltonian reads: 



In TJho the atoms interact with a contact potential modeled by the Dirac (^-function 
with Xij = Xi — Xj denoting the relative coordinate of the i-th and j-th atom. 
The lengths and the energies are scaled by an and TjWii, respectively. Thus, the 
single remaining parameter is the rescaled interaction strength g — n 9 J" a ^ which is 
controlled either by tuning ao via magnetic Feshbach resonances or a± by modifying 
the transversal confinement. This system has been studied for few atoms also in 
[1] [19j [20] . Here we extend the analytical approach we have introduced in [T] to a 
setup involving a larger number of atoms. 

Apart from this single-site prototype system, we are also interested here in finite 
multi-well systems. These can be prepared experimentally e.g. by using more than one 
counter-propagating laser beams (standard optical lattice technique) and forming a 




^ 9 2 1 * 



i=l i t=l 
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superlattice of copies of small finite lattices [5] , or from the beam waist which always 
produces an approximate harmonic confinement resulting in a concentration of the 
cloud density in the few central wells of the potential. The standard optical lattice 
shows a potential profile of the form Vq sin 2 (^p± + cf\ with a superimposed harmonic 
confinement. The rescaled ID Hamiltonian reads: 

1 N rfi N 

ffte= iEp + E 7 ° sin2 + *) + 9 £ s 

i— 1 4 2—1 2<^' 

where the lengths are scaled by the lattice constant a and the energies by the 

ft 2 

recoil energy Er = We further confine this infinite lattice system to a finite 

region L G [-5/275/2] (with <p = ir/2) and L € [-7/2,7/2] (with = 0) 
imposing hard-wall boundary conditions at the edges, such that two and three wells 
are isolated, respectively. This way we encounter simple models of finite multi-well 
systems, a double and a triple well potential. The Hamiltonian is characterized by 
two parameters: the rescaled interaction strength g — and the potential depth 
Vo (in units of Er). We will consider here rather deep lattices (Vq = 40) such that 
the on-site effects can be demonstrated and also a tight-binding approximation with 
well- localized Gaussian funtions at each well, which we will introduce next, can be 
considered as a good approach. 



2.2. Ansatz 

In [I] we have proposed a correlated pair wave function (CPWF) inspired from the idea 
that the analytical solution for a single pair in the harmonic trap can be extended to 
the many-body system by taking products of exact two-body functions. The contact 
interaction may be then adequately addressed, if the discontinuity that it causes is 
imposed on each pair of atoms in the ensemble, in a similar way as for a single pair. 
The CPWF reproduces the two exactly solvable limits of zero and infinite interaction 
strength (Tonks- Girardeau gas) for an arbitrary number of atoms in the harmonic 
trap. We have already indicated in [I] that the CPWF can be used for Monte Carlo 
simulations since it has the form of a Bijl-Jastrow function, but with a different 
approach for the two-body part than the typically used trigonometric functions |15) . 
In this paper we further develop this idea applying it also to higher number of particles, 
and we extend it to the case of a finite-size multi-well potential. 

For a VMC calculation it is crucial to have a good functional form of the trial 
function which can then be used also as a guiding function for DMC simulations 
driving the walkers to positions where, for physical reasons, it is more possible for the 
particles to be placed. For the implementation of DMC we follow similar prescriptions 
like those explained in |15j . The typical trial-guiding function used in these approaches 
use are of the form of Bijl-Jastrow functions, with a single particle part (SPP) which 
accounts for the form of the external potential and an interaction part (IP) which 
accounts for the collision. For a ID harmonic trap the standard form reads: 



^=ri e ^ /2 n cos 



k 3 



(i) 



^>SPP Tp IP 



The form of the SPP corresponds obviously to the ground state of a single particle 
in a harmonic trap (setting also the presumably variational parameter /3 = 1). For the 
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IP the standard Bijl-Jastrow type of functions are used, and an interesting remark is 
that the form cos \k (\xij | — 4)1 is actually the solution of the two-particle Lieb-Liniger 
gas [2] (in homogeneous space) of length L with k determined from the boundary 
condition (or Bethe condition as it is often called) at the contact point X{j = 
[discontinuity of the derivative: 2^(0) = gipij(0) =>■ g = fctan(^)]. For long 
range interparticle distances > L one takes a constant value for the IP. Notice 
that this approach contains one variational parameter L (if we consider /3 fixed) which 
should be optimized. 

Seen from this perspective, our approach in [TJ implies the replacement of the 

IP with the hypergeometric functions U (— fj 5,-^1 [21], which apart from an 

exponential term (which combines with the center of mass motion), are the analytical 

solution of the relative motion of two particles in a trap {^—-^r + -f- — eyj V* = 

0, where = v + \. The parameter v is a generalized quantum number which 
can take real values (and not only integer as in the case of a harmonic trap without 
interaction) and is determined according to the boundary condition at the collision 
point Xij = 0: 

2tr(i^) 
9 r(^) 

Our Ansatz for the wave-function in the harmonic trap then reads: 

^T = U^ X ' /2 U U (-^ X f) ( 2 ) 

i i<j \ / 

and has no variational parameter in the IP, while the one in the SPP can be fixed 
to [$ = 1 according to the ground state of the Hamiltonian Hf lo with g = (being 
simultaneously correct for the Tonks- Girardeau gas limit g — > oo, see [TJ[22])- This 
form Eq. [2] is equivalent to the one appearing in pQ since the relation D u (x) — 

2%e-^U(-%, \, 4) holds for x > |3T]. 

A remark on the two approaches is in order here: the functional form of the 
standard approach [Eq. ([TJ] can come very close to the one proposed here [Eq. (J3J)], 
since it contains a variational parameter L which when optimized leads to a functional 
form for the IP close to the hypergeometric functions. The advantage of the proposed 
approach is that it can be used without any free parameter. Still one can insert 
one variational parameter in the last argument of the hypergeometric function and 
use it variationally too. Nevertheless, our aim in this work is rather to propose a 
function that can describe correctly the behaviour without free parameters. Since both 
approaches are employing in some sense the two-particle solution in the continuum 
[Eq. ([T])] and in the trap [Eq. ((5J)], they offer a physical picture which generalizes the 
two-body problem to the many-body one. 

A different SPP part is needed for the case of a deep finite multi-well potential. In 
the case of a perturbative lattice one could take the same approach or probably slightly 
[Eq. ([2j)] modified by an inverse lattice term for the SPP like 1 — 7sin 2 (7TXi/2 + <p). 
But for a deep lattice Vo > 4(Er) a radical change of the SPP should be employed 
in terms of the tight-binding model which is valid for this case, since next to next- 
neighbor tunneling is at least one order of magnitude smaller than nearest-neighbor 
tunneling. In a general sense one can write an expansion of Gaussian functions (as an 
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approximation to the localized Wannier functions) localized at the position xj of each 
well j: 



The parameter ft is related to the effective confinement frequency of each well 
(considered as a harmonic oscillator), and in our case can be set to ft = y/Va. The 
coefficients fj arc actually dependent on the interaction between the atoms, and in 
the general case can be derived from a Bose-Hubbard model calculation. For weak 
interaction strengths the atoms tend to be more in the center of the potential (for 
harmonic or hard- wall confinement) while increasing the interaction strength leads us 
to the Mott-insulator phase where the particles are essentially equally populating all 
wells (or forming wedding-cake structures of several Mott shells) [TB]. In the case of 
an optical lattice plus a parabolic confinement the tight-binding model for the single- 
particle problem has been solved in [23] by means of Mathieu functions which offer the 
coefficients for an expansion in terms of localized functions. We note here that this 
can be also considered as a parameter free Ansatz for the SPP provided that we use 
the form in Eq. ([3j with a fixed ft = ypfo [21] ■ Yet this will not cover the full range 
from weak to strong interactions since the distribution of weights will change with 
increasing interaction strength. We take a detour from this by focusing on interaction 
effects that appear on-site like those shown in |17] and not on the redistribution of 
the particles from the superfluid to the Mott insulator state. We will examine here 
cases with equal population in each lattice site or obtain this equal distribution for 
a perturbatively weak interaction strength. We will therefore use equal coefficients 
(fj = 1 for all j) in the expansion Eq. ([3J constructing thus a direct generalization of 
our Ansatz for lattices (for periodic lattices this choice is actually the only reasonable). 
In the next section we will examine the adequacy of such an Ansatz to capture the 
properties of the energy and other observables. 

3. Results 

3.1. Harmonic Trap 

In [1] we have shown that for the harmonic trap the energy calculated from the CPWF 
(Eq. [2j is in good agreement with the one obtained from numerical calculations 
for a few atoms. Exact diagonalisation and Multi-Configurational Time-Dependent 
Hartree, both reliable for a few degrees of freedom, have been used there for the 
comparisons. In this work we employ VMC and DMC methods to calculate the energy 
for more atoms up to the standard occupation of each ID tube in experimental setups 
varying from 10-50 atoms. The DMC method [To] , is reliable for extracting the energy 
of a many-body system. 

In Fig. [1] we show the results for the energy of the ground state as a function 
of the interaction strength g for several numbers of particles in the trap. Shown 
is actually the ground state energy of the interacting many-body system divided 
by the energy corresponding to the Tonks- Girardeau limit (which is equal to that 
of the corresponding system of identical fermions Etg — N 2 /2 [5] [55] [T]) as a 
function of the inverse of the interaction strength 1/g (which in fact is proportional 
to the effective ID scattering length ai£> cx — 1/g [II])- The uppermost curve and 
points for N — 5 atoms (as well as for less atoms-not shown) obtained here from 
the DMC calculation and numerical integrations using the functional form of the 




(3) 
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Figure 1. Ground state energies (E) divided by the corresponding value 
of the Tonks-Girardeau limit (Etg) m the harmonic trap as a function of the 
inverse interaction strength (g = g^jj /"hujnau ). Several cases with respect to the 
number of atoms are shown, comparing results obtained by DMC and numerical 
integration of the analytical formula CPWF. 



CPWF Ansatz respectively, confirm the statements in [I] that the CPWF describes 
accurately not only qualitatively but also quantitatively the crossover from weak to 
strong interactions. The agreement is very good allover this crossover, and very 
accurate especially for very weak g — > and very strong interactions close to the 
resonance (TG limit) g — > oo. As we have already shown in [I] the correlated-pair 
wave function represents the exact ground state of the system for the non-interacting 
and TG limit for arbitrary number of atoms. 

For N = 10 or less (see Fig. [lj our Ansatz is still accurate up to 2% error, and 
can be even improved (even below 1%) by VMC letting the parameter j3 of Eq. [5] 
being optimized (this usually leads to lower values of /3) . For an increasing number of 
atoms (compare curves and points in Fig. Q] for N = 20, 30, 50) the intermediate 
interaction regime shows deviations of the Ansatz from the DMC energy curves. 
Still the agreement remains good for strong and weak interaction strengths, and the 
qualitative behaviour throughout the crossover from non-interacting to fermionization 
is captured from the Ansatz. 

3.2. Double and triple well 

The generalization of the CPWF to multi-well systems using the corresponding single- 
particle part (Eq. [3J, is also compared in the following with DMC results. We consider 
here cases that are also relevant for cold-atom experiments, i.e., two to four particles 
per well. 

In Fig. [2] the curves obtained from the Ansatz and the DMC are presented with 
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Figure 2. Ground state energies (E) divided by the corresponding value of the 
Tonks-Girardeau limit (-Etg) f° r ( a ) a double and (b) a triple well as a function of 
the inverse interaction strength (g = gm/Ejia where En is the recoil energy and 
a the lattice constant). Several cases with respect to the number of atoms per well 
(pw) are shown, comparing results obtained by DMC and numerical integration 
of the analytical formula CPWF. 
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Figure 3. One-body density functions p(x) with the lattice constant a as length 
scale for (a) four and (b) six bosons in the double well for several interaction 
strengths covering the crossover from weak to strong couplings. 

respect to the same variables as in Fig. Q] for the double [Fig. [2] (a)] and the triple 
[Fig. [2] (b)] well potentials. The agreement between energies obtained by numericallly 
intergrating the CPWF and DMC calculations is very good for the case of two atoms 
per site both for the double and the triple well potentials. This is to be expected since 
the correlated-pair wave function is actually the exact one in the case of a harmonic 
trap for two atoms [TH [T] . However for three and four atoms per well the deviations 
start to be substantial at intermediate interaction strength (g m 1) and become even 
stronger as we approach the resonance (see Fig. [2]). The reason for these deviations 
lies in the nature of multi-well potentials. As long as the atoms stay in low energy 
states, which is the case for weak interactions the harmonic approach of each well 
remains a good approximation. However, as the interaction strength increases and 
the atoms come energetically closer to the threshold of the barriers of the potential, 
each well becomes more and more anharmonic, as well as the tunneling coupling to 
the neighbouring wells increases. This brcakes the validity of a harmonic approach 
and produces deviations of the Ansatz which still though works qualitatively and as 
we will see next captures important properties of the on-site behaviour beyond the 
weak coupling and the Bose- Hubbard regime. Another important difference compared 
to the case of the harmonic trap is that here the fermionization limit {g — > oo) is 
not reproduced exactly by the Ansatz, due to the anharmonicity explained above. 
Therefore one of the major advantages that our Ansatz proposed in [1] possesses in the 
harmonic trap, namely reproducing the exact behaviour in the two limiting cases (non- 
and infinitely interacting ensemble) , is missing here. Yet the predictions of the Ansatz 
are slightly better if one treats the parameter j3 variationally (VMC calculations show 
a correction up to 1 — 2% in the error), but certainly there is space for improvement 
combining probably Bose-Hubbard and exact solutions for the particular profile of the 
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Figure 4. Two-body density function with the lattice constant a as length scale 
for 4 atoms in the double well for (a) g=2.0 (b) g=5.0 (c) g=10.0 (d) g=20.0 
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multi-well potential in the analytical description. 
3.3. On-site effects 

In general as shown in several publications [501 02] j m ID the density distribution 
of atoms according to the Bose-Hubbard model, does not capture the plethora 
of phenomena occuring beyond standard Mott-insulator. Especially for strong 
interactions interesting on-site effects appear in cases where more than one atom are 
localized in the same well. The density there may broaden or even acquire two or 
more peaks due to the repulsion of the atoms on-site (see refs. |17j). Our Ansatz here 
captures these effects due to the two-body part with the correlated-pair wave function 

The one-body density function p(x) — f_ ... f_ \^(x, X2, x^)\ 2 dx2-..dxN 
(with ty(xi,X2, ...,xn) normalized so that also / p(x)dx = 1) for the cases (a) four 
and (b) six atoms in the double well is shown in Fig. [3] The density funtion at 
each well becomes broader (see g = 2.0, 5.0) acquires a small plateau (g = 10.0) 
and finally develops a number of maxima corresponding to the number of atoms 
localized per site. Similar behaviour has been observed e.g. in [TTJ [50] being 
qualitatively in agreement with our Ansatz. In Fig. @] the two-body density function 
p2{xi,X2) = ■•■ J_ |^| dxs-.-dxN or pair-correlation function of 4 atoms in the 
double well is plotted. For weak interactions [g — 2.0 Fig. 2] (a)] we have only 4 distinct 
peaks due to the double well potential. But as the interaction strength increases 
[g = 5.0 Fig. [4] (b)] the diagonal x\ = X2 tends to deplete (the atoms avoid being 
at the collision points). For stronger interactions [g = 10.0 Fig. [4] (c)] even the off- 
diagonal peaks broaden while close to fermionization [g = 20.0 Fig. [4] (d)] there are 
additional peaks appearing at the off-diagonal spots. Similar effects where shown 
in [201 117j and are in qualitative agreement with those obtained by our analytical 
function. 

4. Conclusions and outlook 

We have studied the crossover from weak to strong interactions of bosonic systems 
in a harmonic trap and in finite multi-well potentials comparing two approaches: 
the numerically exact Quantum Monte Carlo method and the correlated-pair wave 
function approach. Our Ansatz for the harmonic trap introduced and studied in [1] 
for few-body ensembles, is here tested for much larger systems by comparing the 
predictions for the energy with Diffusion Monte Carlo calculations. A good agreement 
all over the interaction regime but particularly close to the Tonks-Girardeau gas 
and for very weak interactions is demonstrated. For larger ensembles there exist 
deviations for the intermediate interaction regime. For finite multi-well systems, 
namely double and triple wells, we have introduced here a modification of the single- 
particle part of the Ansatz in terms of localized functions. We have shown for the 
energy (different loading of the potential from two to four atoms per well) that there 
is qualitative agreement but also strong quantitative deviations as we approach the 
resonance (g — > oo) since the harmonic-approach for each well becomes invalid for 
energetically excited atoms. Still important effects for strong interactions on the one- 
and two-body density functions are captured especially concerning on-site features 
like broadening and appearance of additional peaks. The extension of this idea to 
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fermions or mixtures [2 5) and the improvement of it to describe more accurately the 
full crossover may represent a valuable input to relevant experimental studies. 

Acknowledgements The authors are thankful to Michael Haas and Riidiger Schmitz 
for their help on the implementation of the Quantum Monte Carlo code. 

[1] I. Brouzos and P. Schmelcher, Phys. Rev. Lett. 108, 045301 (2012). 
[2] E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 

[3] Y. Hao, Y. Zhang, J. Q. Liang, and S. Chen, Phys. Rev. A 73, 063617 (2006). 
[4] M. Girardeau, J. Math. Phys. 1, 516 (1960). 

[5] T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004); B. Paredes, A. Wideral, 
V. Murg, O. Mandel, S. Foiling, I. Cirac, G. V. Shlyapnikov, T. W. Hansen, and I. Bloch, 
Nature 429, 277 (2004). 

[6] E. Haller, M. Gustavsson, M. J. Mark, J. G. Danzl, R. Hart, G. Pupillo, H. C. Nagerl, Science 
325, 1224 (2009). 

[7] E. Haller, R. Hart, M.J. Mark, J.G. Danzl, L. Reichsllner, M. Gustavsson, M. Dalmonte, G. 
Pupillo, H.-C. Nagerl, Nature 466, 597 (2010); H. P. Biichler, G. Blatter, and W. Zwerger, 
Phys. Rev. Lett. 90, 130401 (2003) 
[8] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwck, J. Eisert and I. Bloch, Nature 

Physics 8, 325 (2012). 
[9] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008). 
[10] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga , Rev. Mod. Phys. 82, 1225 (2010). 
[11] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). 

[12] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. Gillen, S. Foelling, L. Pollet, M. Greiner, 
Science 329, 547 (2010). 

[13] F. Serwane, G. Ziirn, T. Lompe, T. B. Ottenstein, A. N. Wenz and S. Jochim Science 332, 6027 
(2011). 

[14] G. Ziirn, F. Serwane, T. Lompe, A. N. Wenz, M. G. Ries, J. E. Bohn and S. Jochim, Phys. Rev. 

Lett. 108, 075303 (2012). 
[15] D. Blume and Chris H. Greene, Phys. Rev. A 63, 063601 (2001); D. Blume Phys. Rev. A 66, 

053613 (2002); G. E. Astrakharchik and S. Giorgini, Phys. Rev. A 66, 053614 (2002); G. E. 

Astrakharchik and S. Giorgini, Phys. Rev. A 68, 031602(R) (2003); G. E. Astrakharchik, 

D. Blume, S. Giorgini, and B. E. Granger Phys. Rev. Lett. 92, 030402 (2004); G. E. 
Astrakharchik, D. Blume, S. Giorgini and B. E. Granger, J. Phys. B 37 , 205 (2004); G. 

E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini Phys. Rev. Lett. 95, 190407 
(2005). 

[16] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989); 
D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 
(1998); M. Greiner, O. Mandel, T. Esslinger, T.W. Hansen, and I. Bloch, Nature 415, 39 
(2002); 

[17] I. Brouzos, S. Zollner, and P. Schmelcher, Phys. Rev. A 81, 053613 (2010); O. E. Alon, A. 

I. Streltsov, and L. S. Cederbaum Phys. Rev. Lett. 95, 030405 (2005); D. S. Liihmann, K. 

Bongs, K. Sengstock, and D. Pfannkuche, Phys. Rev. A 77, 023620 (2008), S. Will, T. Best, 

U. Schneider, L. Hackermller, D.S. Liihmann and I. Bloch, Nature 465, 197 (2010) 
[18] T. Busch, B. G. Englert, K. Rzaewski, and M. Wilkens, Found. Phys. 28, 549 (1998). 
[19] F. Deuretzbacher, K. Bongs, K. Sengstock, and D. Pfannkuche, Phys. Rev. A 75, 013614 (2007). 
[20] S. Zollner, H.-D. Meyer, and P. Schmelcher, Phys. Rev. A 74, 063611, 053612 (2006). 
[21] M. Abramowitz and I. A. Stegun Handbook of Mathematical Functions (Dover, New York, 1972). 
[22] M. D. Girardeau, E. M. Wright, and J. M. Triscari, Phys. Rev. A 63, 033601 (2001). 
[23] A. M. Rey, G. Pupillo, C. W. Clark, and C. J. Williams. Phys. Rev. A 72, 033616 (2005). 
[24] According to 1231 fj are the Fourier coefficients of the Mathieu function ce(x,—q): fj = 

~ J ce (x, —q) cos (2jx) dx. The parameter q = where: 

J = - J w( Xj ) \-]-j^ + V A^{wx/2)\ u^-i) (4) 

characterizes the system, where w(x) are the localized Wannier functions for each well . If 
Gaussian profile functions are used for w(x) like in Eq.|3]one obtains 



q(V ) = -VS^V 1/4 e- 2 V^o 
[25] I. Brouzos and P. Schmelcher, larXiv: 1209.28911 



1 + e 8 vT¥ y/v 



(5) 



